#Ejercicios {.tabset .tabset-fade .tabset-pills}

Capitulo IV

###Punto 10

(a) Produzaca algunos resumenes numéricos y gráficos de los datos Weekly.¿ parace que hay patrones?.

Resumen de variables:

library(ISLR)
data("Weekly")
summary(Weekly)
##       Year           Lag1               Lag2               Lag3         
##  Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
##  1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580  
##  Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410  
##  Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472  
##  3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
##  Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
##       Lag4               Lag5              Volume       
##  Min.   :-18.1950   Min.   :-18.1950   Min.   :0.08747  
##  1st Qu.: -1.1580   1st Qu.: -1.1660   1st Qu.:0.33202  
##  Median :  0.2380   Median :  0.2340   Median :1.00268  
##  Mean   :  0.1458   Mean   :  0.1399   Mean   :1.57462  
##  3rd Qu.:  1.4090   3rd Qu.:  1.4050   3rd Qu.:2.05373  
##  Max.   : 12.0260   Max.   : 12.0260   Max.   :9.32821  
##      Today          Direction 
##  Min.   :-18.1950   Down:484  
##  1st Qu.: -1.1540   Up  :605  
##  Median :  0.2410             
##  Mean   :  0.1499             
##  3rd Qu.:  1.4050             
##  Max.   : 12.0260

Matriz de gráfico de dispersión:

pairs(Weekly)

En la primera matriz muestra un resumen de las variables a trabajar, observamos que todas estan a una misma escala, ecepto la varible dirección que cuenta con dos niveles (Alto y bajo). En la segunda matriz las varibles lag1,lag2,lag3,lag4,lag5 y volumen, son variables que entre ellas aparentemente las observaciones se concentrar en el centro; podemos notar que las varible año y variable volume tiene un comportamiento de un polinomio de grado dos acsendente mientras pasa el tiempo.

Matriz de correlaciones

cor(Weekly[,-9])
##               Year         Lag1        Lag2        Lag3         Lag4
## Year    1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923
## Lag1   -0.03228927  1.000000000 -0.07485305  0.05863568 -0.071273876
## Lag2   -0.03339001 -0.074853051  1.00000000 -0.07572091  0.058381535
## Lag3   -0.03000649  0.058635682 -0.07572091  1.00000000 -0.075395865
## Lag4   -0.03112792 -0.071273876  0.05838153 -0.07539587  1.000000000
## Lag5   -0.03051910 -0.008183096 -0.07249948  0.06065717 -0.075675027
## Volume  0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617
## Today  -0.03245989 -0.075031842  0.05916672 -0.07124364 -0.007825873
##                Lag5      Volume        Today
## Year   -0.030519101  0.84194162 -0.032459894
## Lag1   -0.008183096 -0.06495131 -0.075031842
## Lag2   -0.072499482 -0.08551314  0.059166717
## Lag3    0.060657175 -0.06928771 -0.071243639
## Lag4   -0.075675027 -0.06107462 -0.007825873
## Lag5    1.000000000 -0.05851741  0.011012698
## Volume -0.058517414  1.00000000 -0.033077783
## Today   0.011012698 -0.03307778  1.000000000

(b) Utilice el conjunto de datos completo para realizar una regresion logística con Direction como la respuesta y la cinco varibles de retraso mas Volume como predictoras. Use la funcion summary para imprimir los resultados. ¿algunos de los predictores parece ser estadisticamente significativos? si es asi ¿cuales son estas variables?

Resumen Regresión logística

mod1 <- glm( Direction ~ Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data = Weekly,family = "binomial")
summary(mod1)
## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = "binomial", data = Weekly)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6949  -1.2565   0.9913   1.0849   1.4579  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.26686    0.08593   3.106   0.0019 **
## Lag1        -0.04127    0.02641  -1.563   0.1181   
## Lag2         0.05844    0.02686   2.175   0.0296 * 
## Lag3        -0.01606    0.02666  -0.602   0.5469   
## Lag4        -0.02779    0.02646  -1.050   0.2937   
## Lag5        -0.01447    0.02638  -0.549   0.5833   
## Volume      -0.02274    0.03690  -0.616   0.5377   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1496.2  on 1088  degrees of freedom
## Residual deviance: 1486.4  on 1082  degrees of freedom
## AIC: 1500.4
## 
## Number of Fisher Scoring iterations: 4

Podemos observar que la variable Lag2 con un P-valor de 0.0296 es estadisiticamente significativa con un nivel de confianza de 0.05.

(c) Compute la matrix de confusión y la fracción general de las predicciones correctas. Explica qué te dice la matriz de confusión logística.

Tabla de confusión

prediction <- mod1$fitted.values
pred<- rep("Dow",length(prediction))
pred[prediction > 0.5]<- "Up"
table(pred,Weekly$Direction)
##      
## pred  Down  Up
##   Dow   54  48
##   Up   430 557

Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{54+557}{1089} = 0.5610 \times 100\) = 56.10% del tiempo . La tasa de error con los datos de entrenamiento es del 100%-56.10% = 43.8%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{557}{48+557}=0.92\times100%\)=92% de la veces, tambien cuando el mercado baja el modelo acierta un \(\frac{54}{54+430}=0.1115\times 100\) = 11.15% del tiempo.

(d) Ahora ajuste el modelo de regresión logística utilizando un período de datos de entrenamiento de 1990 a 2008, con Lag2 como el unico predictor. Calcule la matriz de confusión y la fracción general de las predicciones correctas para los datos retenidos (Es decir, los datos de 2009 y 2010)

Resumen Modelo con lag2.

train <- subset.data.frame(x = Weekly,subset = Year < 2009)
test2009_2010 <- subset.data.frame(Weekly,subset = Year >=2009)
mod2 <- glm(Direction ~ Lag2 ,data = train,family = binomial)
summary(mod2)
## 
## Call:
## glm(formula = Direction ~ Lag2, family = binomial, data = train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.536  -1.264   1.021   1.091   1.368  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.20326    0.06428   3.162  0.00157 **
## Lag2         0.05810    0.02870   2.024  0.04298 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1354.7  on 984  degrees of freedom
## Residual deviance: 1350.5  on 983  degrees of freedom
## AIC: 1354.5
## 
## Number of Fisher Scoring iterations: 4

Tabla de confusión

prediction2 <- predict(object = mod2,newdata = test2009_2010,type = "response")
pred2<- rep("Dow",length(prediction2))
pred2[prediction2 > 0.5]<- "Up"
table(pred2,test2009_2010$Direction)
##      
## pred2 Down Up
##   Dow    9  5
##   Up    34 56

Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{9+56}{104}=0.625\times 100\)= 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-62.5% = 37.5%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{56}{5+56}=0.9180\times100%\)=91.8% del tiempo, tambien cuando el mercado baja el modelo acierta un \(\frac{34}{34+9}=0.79\times 100\) = 79.6% del tiempo.

(e) Repita (d) usando LDA

Resumen modelo LDA:

library(MASS)
modlda <- lda(Direction~Lag2,data = train)
modlda
## Call:
## lda(Direction ~ Lag2, data = train)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581
## 
## Coefficients of linear discriminants:
##            LD1
## Lag2 0.4414162

Tabla de confusión:

prediction_LDA <- predict(object = modlda,newdata = test2009_2010)
table(prediction_LDA$class,test2009_2010$Direction)
##       
##        Down Up
##   Down    9  5
##   Up     34 56

Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{9+56}{104}=0.625\times 100\)= 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-62.5% = 37.5%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{56}{5+56}=0.9180\times100%\)=91.8% de la veces, tambien cuando el mercado baja el modelo acierta un \(\frac{34}{34+9}=0.79\times 100\) = 79.6% del tiempo.

(f) Repita (d) usando QDA

Resumen modelo QDA

mod_QDA <-  qda(Direction~Lag2,data = train)
mod_QDA
## Call:
## qda(Direction ~ Lag2, data = train)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581
prediction_QDA <- predict(mod_QDA,newdata = test2009_2010)
table(prediction_QDA$class,test2009_2010$Direction)
##       
##        Down Up
##   Down    0  0
##   Up     43 61

Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{0+61}{104}=0.5865 \times 100\)= 58.65% del tiempo. La tasa de error con los datos de prueba es del 100%-58.65% = 41.35%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{61}{61}=1\) osea el 100% de la veces, Se resalta este modelo de anialisis de discriminación cudratica elige una dirección alta todo el tiempo.

(g) Repita (d) usando Knn con k=1.

Table de confusión:

library(class)
train.x <- as.matrix(train$Lag2)
test.x <- as.matrix(test2009_2010$Lag2)
train.Direction <- train$Direction
set.seed(1)
prediction_knn <- knn(train = train.x,test = test.x,cl = train.Direction,k = 1)
table(prediction_knn,test2009_2010$Direction)
##               
## prediction_knn Down Up
##           Down   21 30
##           Up     22 31

Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-56.10% = 43.9%. Tambien apreciamos que cuando el mercado sube, El modelo acierta 50.81% de la veces, tambien cuando el mercado baja el modelo acierta un 48.83% del tiempo.

(h) ¿ cuál de estos métodos parece proporcionar los mejores resultados con estos datos ?

Con los modelos anteriores, vemos que el modelo de regresión logístico y LDA, La tasa de error es mínima. tambien para los modelos QDA y KNN un poco menores.

(i) Experimente con diferentes combinaciones de predictores, incluyendo posibles transformaciones, para cada uno de los métodos, reporte la variables,método y la matriz de confusión asociada que parece proporcionar los mejores resultados en los mejores resultados en los datos retenidos.Tenga encuenta también debe experimentar con los valores para K en la clasificación con KNN.

La variable mas significativa es Lag2 como se mostro en los anteriores puntos para tener una interacción con las segunda varible significativa.

library(MASS)
library(class)

mod_inter_1 <- glm(Direction~Lag2:Lag1,family = binomial,data = train)
prediction_inter_1 <- predict(mod_inter_1,test2009_2010,type = "response")
pred_inter_1 <- rep("Dow",length(prediction_inter_1))
pred_inter_1[prediction_inter_1 > 0.5] <- "Up"


mod_lda_inter <-lda(Direction~Lag2:Lag1,data = train)

prediction_LDA_inter <- predict(mod_lda_inter,test2009_2010)


mod_QDA_inter <- qda(Direction ~ Lag2+sqrt(abs(Lag2)),train)
pred_QDA_inter <- predict(mod_QDA_inter,test2009_2010)

set.seed(0511)
#knn k=10
prediction_knn1 <- knn(train = train.x,test = test.x,cl = train.Direction,k = 10)
#knn=100
prediction_knn2 <- knn(train = train.x,test = test.x,cl = train.Direction,k = 100)

Tablas de modelos realizados con interacciones

Modelo predicción correcta tasa de presición cuando el mercado aumenta
regresión logistica con interacción 57.69% 98.36%
LDA con interacción 57.69% 98.36%%
QDA con \(\sqrt(abs(Lag2))\) 57.69% 78.68%
knn con K=10 57.69% 68.85%
Knn con k=100 56.73% 80.32%

De esta tabla se puede concluir que el modelo de regresión logística y LDA tienen el mejor rendimiento en tasas de error de prueba.

Punto 11

  • En este problema, desarrollará un modelo para predecir si un automóvil determinado obtiene un consumo de combustible alto o bajo en función del conjunto de datos automático.

(a) Cree una variable binaria, mpg01, que contenga un 1 si mpg contiene un valor por encima de su mediana, y un 0 si mpg contiene un valor por debajo de su mediana. Puede calcular la mediana utilizando la función median(). Tenga en cuenta que puede resultarle útil utilizar la función data.frame() para crear un único conjunto de datos que contenga tanto mpg01 como las otras variables automáticas.

library(ISLR)
data(Auto)
n <- length(Auto$mpg)

mpg01 <- ifelse( Auto$mpg > median(Auto$mpg),1,0)
Auto <- data.frame(Auto,mpg01)
Auto$mpg01 <- as.factor(Auto$mpg01)
knitr::kable(head(Auto),caption = "Cabeza de base de datos Auto con variable mpg01")
Cabeza de base de datos Auto con variable mpg01
mpg cylinders displacement horsepower weight acceleration year origin name mpg01
18 8 307 130 3504 12.0 70 1 chevrolet chevelle malibu 0
15 8 350 165 3693 11.5 70 1 buick skylark 320 0
18 8 318 150 3436 11.0 70 1 plymouth satellite 0
16 8 304 150 3433 12.0 70 1 amc rebel sst 0
17 8 302 140 3449 10.5 70 1 ford torino 0
15 8 429 198 4341 10.0 70 1 ford galaxie 500 0

(b) Explore los datos gráficamente para investigar la asociación entre mpg01 y las otras características. ¿Cuál de las otras características parece más útil para predecir mpg01? Los diagramas de dispersión y los diagramas de caja pueden ser herramientas útiles para responder a esta pregunta. Describe tus hallazgos.

Correlaciones entre variables:

Auto1 <- Auto
Auto1$mpg01 <- as.numeric(Auto1$mpg01)
cor(Auto1[,-9])
##                     mpg  cylinders displacement horsepower     weight
## mpg           1.0000000 -0.7776175   -0.8051269 -0.7784268 -0.8322442
## cylinders    -0.7776175  1.0000000    0.9508233  0.8429834  0.8975273
## displacement -0.8051269  0.9508233    1.0000000  0.8972570  0.9329944
## horsepower   -0.7784268  0.8429834    0.8972570  1.0000000  0.8645377
## weight       -0.8322442  0.8975273    0.9329944  0.8645377  1.0000000
## acceleration  0.4233285 -0.5046834   -0.5438005 -0.6891955 -0.4168392
## year          0.5805410 -0.3456474   -0.3698552 -0.4163615 -0.3091199
## origin        0.5652088 -0.5689316   -0.6145351 -0.4551715 -0.5850054
## mpg01         0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566
##              acceleration       year     origin      mpg01
## mpg             0.4233285  0.5805410  0.5652088  0.8369392
## cylinders      -0.5046834 -0.3456474 -0.5689316 -0.7591939
## displacement   -0.5438005 -0.3698552 -0.6145351 -0.7534766
## horsepower     -0.6891955 -0.4163615 -0.4551715 -0.6670526
## weight         -0.4168392 -0.3091199 -0.5850054 -0.7577566
## acceleration    1.0000000  0.2903161  0.2127458  0.3468215
## year            0.2903161  1.0000000  0.1815277  0.4299042
## origin          0.2127458  0.1815277  1.0000000  0.5136984
## mpg01           0.3468215  0.4299042  0.5136984  1.0000000

Grafico de dispersion

pairs(Auto)

attach(Auto)
## The following object is masked _by_ .GlobalEnv:
## 
##     mpg01
layout(rbind(c(1,1,1,2,2,2,3,3,3),c(4,4,4,5,5,5,6,6,6)))
boxplot(cylinders~mpg01,data = Auto,main="Cilindros vs mpg01")
boxplot(displacement ~mpg01,data = Auto,main="Desplazamiento vs mpg01")
boxplot(horsepower ~mpg01,data = Auto,main="Caballos de fuerza vs mpg01")
boxplot(weight ~mpg01,data = Auto,main="Peso del vehiculo vs mpg01")
boxplot(acceleration ~mpg01,data = Auto,main="Aceleración vs mpg01")
boxplot(year ~mpg01,data = Auto,main="Modelo vs mpg01")

(c) Divida los datos en un conjunto de entrenamiento y un conjunto de prueba.

train <- sample(seq(length(Auto$mpg01)),length(Auto$mpg01)*0.70,replace = FALSE)
data_train <- Auto[train,]
data_test <- Auto[-train,]
Entranamiento Prueba Total
274 118 392
————- ——- —–

(d) Realice LDA en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?

require(MASS)
fit.lda <- lda(mpg01 ~ cylinders+displacement+horsepower+weight,Auto)
fit.lda
## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto)
## 
## Prior probabilities of groups:
##   0   1 
## 0.5 0.5 
## 
## Group means:
##   cylinders displacement horsepower   weight
## 0  6.765306     273.1582  130.11224 3620.403
## 1  4.178571     115.6658   78.82653 2334.765
## 
## Coefficients of linear discriminants:
##                        LD1
## cylinders    -0.4708126926
## displacement -0.0014339977
## horsepower    0.0044241441
## weight       -0.0009846242
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$mpg01)
##    
##      0  1
##   0 51  3
##   1 11 53
mean(pred.lda$class != data_test$mpg01)
## [1] 0.1186441

Con los resultados anteriores podemos concluir que la tasa de error es del 11.01%.

(e) Realice QDA en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?

fit.qda <- qda(mpg01~cylinders+displacement+horsepower,data = data_train)
fit.qda
## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower, data = data_train)
## 
## Prior probabilities of groups:
##         0         1 
## 0.4890511 0.5109489 
## 
## Group means:
##   cylinders displacement horsepower
## 0  6.783582     278.5970   132.8881
## 1  4.185714     117.8071    80.1500
pred.qda <- predict(fit.qda,data_test)
table(pred.qda$class,data_test$mpg01)
##    
##      0  1
##   0 53  3
##   1  9 53
mean(pred.qda$class != data_test$mpg01)
## [1] 0.1016949

Con los resultados anteriores podemos concluir que la tasa de error es del 11.01%.

(f) Realice una regresión logística en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?

fit.glm <- glm(mpg01 ~ cylinders+displacement+horsepower+weight+acceleration,family = binomial,data_train)
summary(fit.glm)
## 
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower + 
##     weight + acceleration, family = binomial, data = data_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4925  -0.1566   0.1188   0.3611   3.2309  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  11.730488   3.276787   3.580 0.000344 ***
## cylinders    -0.136576   0.404783  -0.337 0.735811    
## displacement -0.007845   0.009597  -0.817 0.413649    
## horsepower   -0.033194   0.024676  -1.345 0.178556    
## weight       -0.002468   0.001116  -2.212 0.026996 *  
## acceleration  0.038322   0.159330   0.241 0.809930    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 379.71  on 273  degrees of freedom
## Residual deviance: 147.81  on 268  degrees of freedom
## AIC: 159.81
## 
## Number of Fisher Scoring iterations: 7
pr<- predict(fit.glm,data_test,type = "response")
pred.glm <- rep(0,length(pr))
pred.glm[pr>0.5]<-1
table(pred.glm,data_test$mpg01)
##         
## pred.glm  0  1
##        0 52  3
##        1 10 53
mean(pred.glm != data_test$mpg01)
## [1] 0.1101695

Con los resultados anteriores podemos concluir que la tasa de error es del 11.86%.

(g) Realice KNN en los datos de entrenamiento, con varios valores de K, para predecir mpg01. Use solo las variables que parecían más asociadas con mpg01 en (b). ¿Qué errores de prueba obtienes? ¿Qué valor de K parece tener el mejor rendimiento en este conjunto de datos?

require(ggplot2)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'Auto':
## 
##     mpg
require(factoextra)
## Loading required package: factoextra
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
require(class)
require(caret)
## Loading required package: caret
## Loading required package: lattice
set.seed(0511)
ctrl <- trainControl(method="repeatedcv",repeats = 3) 
knnFit <- train(mpg01 ~ cylinders+displacement+horsepower+weight+acceleration, data =data_train[,-c(1,7,8,9)], method = "knn", trControl = ctrl, preProcess = c("center","scale"), tuneLength = 20)
knnFit
## k-Nearest Neighbors 
## 
## 274 samples
##   5 predictor
##   2 classes: '0', '1' 
## 
## Pre-processing: centered (5), scaled (5) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 247, 247, 247, 247, 246, 246, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    5  0.9063492  0.8123034
##    7  0.9003086  0.8002292
##    9  0.8966931  0.7931661
##   11  0.9039683  0.8077058
##   13  0.9039683  0.8077058
##   15  0.9039683  0.8077058
##   17  0.9039683  0.8076792
##   19  0.9039683  0.8076792
##   21  0.9015432  0.8028595
##   23  0.9015432  0.8028595
##   25  0.9015432  0.8028595
##   27  0.9015432  0.8028595
##   29  0.9015432  0.8028595
##   31  0.9015432  0.8028595
##   33  0.9015432  0.8028595
##   35  0.9015432  0.8028595
##   37  0.9015432  0.8028595
##   39  0.9027337  0.8052404
##   41  0.9027337  0.8052404
##   43  0.9015432  0.8028595
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.

se obtuvo un k optimo de 7 procedemos con el modelo.

train.x <- as.matrix(data_train[,-c(1,7,8,9)])
test.x <- as.matrix(data_test[,-c(1,7,8,9)])
train.y <- as.matrix(data_train$mpg01)
test.y <- as.matrix(data_test$mpg01)
set.seed(0511)
pred.knn <- knn(train = train.x,test = test.x,cl = train.y,k = 7)
table(pred.knn,test.y)
##         test.y
## pred.knn  0  1
##        0 52  4
##        1 10 52
mean(pred.knn !=test.y)
## [1] 0.1186441

con los resultados obtenidos con un K=7 optimo la tasa de error es de 13.55%.

Conclución general: Los mejores modelos son LDA y QDA con una tasa de error de 11.01%

Punto 12

  • Este problema implica escribir funciones.

(a) Escriba una función, Power(), que imprima el resultado de elevar 2 a la tercera potencia. En otras palabras, su función debe calcular \(2^3\) e imprimir los resultados.Sugerencia: recuerde que \(x^a\) eleva x a la potencia a. Use la función print() para generar el resultado.

Power <- function() {
    2^3
}

Power()
## [1] 8

(b) Cree una nueva función, Power2(), que le permita pasar dos números, x y a, e imprima el valor de \(x^a\). Puede hacer esto comenzando su función con la línea

power <- function(){2^3}
power()
## [1] 8

(b) Cree una nueva función, Power2(), que le permita pasar dos números, x y a, e imprima el valor de \(x^a\). Puede hacer esto comenzando su función con la línea

Power2 <- function (x,a){x^a}

Power2(3,8)
## [1] 6561
Power2 <- function(x, a) {
    x^a
}

Power2(3, 8)
## [1] 6561

(c) Usando la función Power2() que acaba de escribir, calcule \(10^3\), \(8^17\) y \(131^3\).

Power2(10,3)
## [1] 1000
Power2(8,17)
## [1] 2.2518e+15
Power2(131,3)
## [1] 2248091

(d) Ahora cree una nueva función, Power3(), que realmente devuelve el resultado x^a como un objeto R, en lugar de simplemente imprimirlo en la pantalla. Es decir, si almacena el valor x^a en un objeto llamado resultado dentro de su función, simplemente puede return() este resultado, utilizando la siguiente línea:

power3 <- function(x,a){
  resultado <- x^a
  return(resultado)
}
Power3 <- function(x , a) {
    result <- x^a
    return(result)
}

(e) Ahora, utilizando la función Power3 (), cree un plot de f(x) = \(x^2\). El eje x debe mostrar un rango de enteros del 1 al 10, y el eje y debería mostrar \(x^2\). Rotule los ejes apropiadamente y use un título apropiado para la figura. Considere mostrar el eje x, el eje y o ambos en la escala logarítmica. Puede hacerlo utilizando log = ‘‘ x’’, log = ‘‘ y’’ o log = ‘‘ xy’’ como argumentos de la función plot()

función a graficar:

\[f(x)=x^2\]

x <- 1:10
plot(x,power3(x,2),log = "xy",xlab = "log de x",ylab = "log de x^2",main = "log de x^2 vs log de x")

x <- 1:10
plot(x, Power3(x, 2), log = "xy", xlab = "Log of x", ylab = "Log of x^2", main = "Log of x^2 vs Log of x")

(f) Cree una función, PlotPower (), que le permite crear una gráfica de x contra x ^ a para un a fijo y para un rango de valores de x. Por ejemplo, si llamas

entonces se debe crear un gráfico con un eje x que tome valores 1,2, …, 10 y un eje y que toma los valores \(1^3\),\(2^3\), …, \(10^3\).

plotpower <- function(x,a){
  plot(x,power3(x,a))
}
plotpower(1:10,3)

PlotPower <- function(x, a) {
    plot(x, Power3(x, a))
}

PlotPower(1:10, 3)

Punto 13

Utilizando el conjunto de datos de Boston, ajuste los modelos de clasificación para predecir si un suburbio determinado tiene una tasa de criminalidad superior o inferior a la mediana. Explore los modelos de regresión logística, LDA y KNN utilizando varios subconjuntos de predictores. Describe tus hallazgos.

library(MASS)
data("Boston")
medianacrim <- median(Boston$crim)
crim01 <- rep(0,length(Boston$crim))
crim01[Boston[1]>medianacrim] <- 1
Boston$crim01 <- crim01
train <- sample(seq(length(Boston$crim01)),length(Boston$crim01)*0.70,replace = FALSE)
data_train <- Boston[train,]
data_test <- Boston[-train,]

Regresion logistica:

fit.glm <- glm(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat,family = binomial,data_train)
summary(fit.glm)
## 
## Call:
## glm(formula = crim01 ~ zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat, family = binomial, data = data_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9848  -0.1797  -0.0017   0.0034   3.3317  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -30.570448   7.823691  -3.907 9.33e-05 ***
## zn           -0.066360   0.039966  -1.660 0.096830 .  
## indus        -0.043417   0.053660  -0.809 0.418451    
## chas          0.812084   0.893650   0.909 0.363495    
## nox          46.839608   8.845996   5.295 1.19e-07 ***
## rm            0.280004   0.563847   0.497 0.619474    
## age           0.022823   0.014516   1.572 0.115881    
## dis           0.483227   0.236176   2.046 0.040752 *  
## rad           0.693487   0.184343   3.762 0.000169 ***
## tax          -0.009323   0.003294  -2.830 0.004649 ** 
## ptratio       0.178923   0.121506   1.473 0.140874    
## black        -0.006888   0.005393  -1.277 0.201544    
## lstat        -0.036400   0.057403  -0.634 0.526006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.70  on 353  degrees of freedom
## Residual deviance: 151.08  on 341  degrees of freedom
## AIC: 177.08
## 
## Number of Fisher Scoring iterations: 9
prob.glm <- predict(fit.glm,data_test,type = "response")
pred.glm <- rep(0,length(prob.glm))
pred.glm[prob.glm > 0.5] <- 1
table(pred.glm,data_test$crim01)
##         
## pred.glm  0  1
##        0 65  9
##        1  9 69
mean(pred.glm != data_test$crim01)
## [1] 0.1184211

con los resultados obtenidos la tasa de error es de 13.15%

Con variables predictora con correlacion fuerte

fit.glm1 <- glm(crim01~indus+nox+age+dis+rad+tax+black+lstat,family = binomial,data_train)
summary(fit.glm1)
## 
## Call:
## glm(formula = crim01 ~ indus + nox + age + dis + rad + tax + 
##     black + lstat, family = binomial, data = data_train)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.96371  -0.22572  -0.01948   0.00597   2.91220  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -23.472808   5.213152  -4.503 6.71e-06 ***
## indus        -0.052668   0.054397  -0.968  0.33293    
## nox          45.719887   8.902831   5.135 2.81e-07 ***
## age           0.024557   0.012887   1.906  0.05671 .  
## dis           0.265804   0.187558   1.417  0.15643    
## rad           0.620643   0.146227   4.244 2.19e-05 ***
## tax          -0.009174   0.003001  -3.057  0.00223 ** 
## black        -0.007900   0.005400  -1.463  0.14348    
## lstat        -0.042885   0.040838  -1.050  0.29366    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 490.70  on 353  degrees of freedom
## Residual deviance: 160.24  on 345  degrees of freedom
## AIC: 178.24
## 
## Number of Fisher Scoring iterations: 8
prob.glm1 <- predict(fit.glm1,data_test,type = "response")
pred.glm1 <- rep(0,length(prob.glm1))
pred.glm1[prob.glm1 > 0.5] <- 1
table(pred.glm1,data_test$crim01)
##          
## pred.glm1  0  1
##         0 67 12
##         1  7 66
mean(pred.glm1 != data_test$crim01)
## [1] 0.125

con estos variables predictoras la tasa de error es de 15.13%.

LDA:

fit.lda <- lda(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat , data = data_train)
fit.lda
## Call:
## lda(crim01 ~ zn + indus + chas + nox + rm + age + dis + rad + 
##     tax + ptratio + black + lstat, data = data_train)
## 
## Prior probabilities of groups:
##         0         1 
## 0.5056497 0.4943503 
## 
## Group means:
##          zn     indus       chas       nox       rm      age      dis
## 0 22.460894  7.143631 0.04469274 0.4708587 6.399620 50.44469 5.153022
## 1  1.268571 15.137771 0.09142857 0.6337086 6.170234 85.79200 2.523521
##         rad      tax  ptratio    black     lstat
## 0  4.201117 311.8156 17.92514 388.2874  9.549944
## 1 14.777143 506.8000 18.96800 319.7955 15.626343
## 
## Coefficients of linear discriminants:
##                   LD1
## zn      -0.0034199390
## indus    0.0102997449
## chas     0.0808703228
## nox      7.2261112533
## rm       0.0612382181
## age      0.0167656971
## dis     -0.0027187266
## rad      0.0970467645
## tax     -0.0018560309
## ptratio -0.0064709707
## black   -0.0005065647
## lstat   -0.0255231030
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$crim01)
##    
##      0  1
##   0 69 15
##   1  5 63
mean(pred.lda$class != data_test$crim01)
## [1] 0.1315789

Con los resultados obtenidos con este metodo, la tasa de error es de 17.10% mas que la regresión logistíca.

Ahora con las variables predictoras con correlación fuerte.

fit.lda1 <- lda(crim01~indus+nox+age+dis+rad+tax+black+lstat , data = data_train)
fit.lda1
## Call:
## lda(crim01 ~ indus + nox + age + dis + rad + tax + black + lstat, 
##     data = data_train)
## 
## Prior probabilities of groups:
##         0         1 
## 0.5056497 0.4943503 
## 
## Group means:
##       indus       nox      age      dis       rad      tax    black
## 0  7.143631 0.4708587 50.44469 5.153022  4.201117 311.8156 388.2874
## 1 15.137771 0.6337086 85.79200 2.523521 14.777143 506.8000 319.7955
##       lstat
## 0  9.549944
## 1 15.626343
## 
## Coefficients of linear discriminants:
##                 LD1
## indus  0.0118027470
## nox    7.2267222513
## age    0.0174936477
## dis   -0.0237077126
## rad    0.0998803157
## tax   -0.0021239136
## black -0.0005370087
## lstat -0.0291598636
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$crim01)
##    
##      0  1
##   0 69 15
##   1  5 63
mean(pred.lda$class != data_test$crim01)
## [1] 0.1315789

Con estas covariables se obtiene una tasa de error de 17.10% igual que la anterior pero por pricipio de parcimonia es mejor este modelo con pocas variables

KNN

set.seed(0511)
data_train$crim01 <- as.factor(data_train$crim01)
ctrl <- trainControl(method="repeatedcv",repeats = 3) 
knnFit <- train(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat, data =data_train, method = "knn", trControl = ctrl, preProcess = c("center","scale"), tuneLength = 20)
knnFit
## k-Nearest Neighbors 
## 
## 354 samples
##  12 predictor
##   2 classes: '0', '1' 
## 
## Pre-processing: centered (12), scaled (12) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 318, 319, 319, 319, 318, 319, ... 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    5  0.9057672  0.8113559
##    7  0.8755556  0.7509298
##    9  0.8576984  0.7151367
##   11  0.8641799  0.7281249
##   13  0.8500000  0.6996443
##   15  0.8396296  0.6788571
##   17  0.8302646  0.6600574
##   19  0.8358201  0.6711574
##   21  0.8320635  0.6635468
##   23  0.8284127  0.6562289
##   25  0.8181746  0.6357683
##   27  0.8144444  0.6283471
##   29  0.8152116  0.6296212
##   31  0.8113492  0.6218266
##   33  0.8114286  0.6218624
##   35  0.8096032  0.6181100
##   37  0.8076984  0.6142139
##   39  0.7888889  0.5764674
##   41  0.7907672  0.5801909
##   43  0.7935714  0.5857363
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.

Un K optimo para este metodo de clasificación es 5.

train.x <- as.matrix(data_train[,-c(1,14)])
test.x <- as.matrix(data_test[,-c(1,14)])
train.y <- as.matrix(data_train$crim01)
test.y <- as.matrix(data_test$crim01)
set.seed(0511)
pred.knn <- knn(train = train.x , test = test.x,cl = train.y,k=5)
table(pred.knn,test.y)
##         test.y
## pred.knn  0  1
##        0 66  7
##        1  8 71
mean(pred.knn !=test.y)
## [1] 0.09868421

Este metodo tiene una tasa de error de 7.2% mucho mejor que los anteriores.

Capitulo VIII

###Punto 7

library(MASS)
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
train <- sample(1:nrow(Boston),length(Boston$crim)*0.70,replace = FALSE)
data.train <- Boston[train,-14]
data.test <- Boston[-train,-14]
y.train <- Boston[train,14]
y.test <- Boston[-train,14]
set.seed(0511)
boston.1 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
                         ntree = 500,mtry = ncol(Boston)-1)

boston.2 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
                         ntree = 500,mtry = (ncol(Boston)-1)/2)

boston.3 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
                         ntree = 500,mtry = sqrt(ncol(Boston)-1))

val <- data.frame(arboles = 1:500)
val["MSE1"] <- boston.1$test$mse
val["MSE2"] <- boston.2$test$mse
val["MSE3"] <- boston.3$test$mse

library(ggplot2)
library(reshape2)
xy <- melt(val,id=c("arboles"))

theme_set(theme_bw())
ggplot(xy)+geom_line(aes(x=arboles,y = value,color=variable))+
  scale_color_manual(name=" m = ",labels=c("p","p/2","\u221Ap"),values = c("orange","red","green"))+labs(x = "Número de árboles",
       y = "Error de clasificación") +
  theme(legend.position="top")

Con esta grafica podemos observar que la prueba MSE es muy alta para un solo arbol, a mededida que aunmenta los arboles el error se estabiliza aproximademente cuando hay 150 arboles, La estaibilización se presenta para los tres valores de “m”. El m que menor reprensenta error es m = \(sqrt(p)\).Ojo que al variar la semilla se podra concluir diferente.

###Punto 8

(a) Divida el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba.

library(ISLR)
set.seed(0511)
train <- sample(1:nrow(Carseats),nrow(Carseats)*0.65,replace = FALSE)
carseats.train <- Carseats[train,]
carseats.test <- Carseats[-train,]

(b) Ajuste un árbol de regresión al conjunto de entrenamiento. Trace el árbol e interprete los resultados. ¿Qué prueba MSE obtienes?

library(rpart)
library(rpart.plot)
carseats.tree <- rpart(Sales~.,data = carseats.train)
rpart.plot(carseats.tree,box.palette = "RdBu",nn=TRUE)

pred <- predict(carseats.tree,carseats.test)
mean((pred - carseats.test$Sales)^2)
## [1] 5.463956

con el anterior resultado concluimos que el MSE es aproximadamente de 5.46.

(c) Utilice la validación cruzada para determinar el nivel óptimo de complejidad del árbol. ¿La poda del árbol mejora la prueba MSE?

printcp(carseats.tree)
## 
## Regression tree:
## rpart(formula = Sales ~ ., data = carseats.train)
## 
## Variables actually used in tree construction:
## [1] Advertising Age         CompPrice   Income      Price       ShelveLoc  
## [7] US         
## 
## Root node error: 2084.7/260 = 8.0181
## 
## n= 260 
## 
##          CP nsplit rel error  xerror     xstd
## 1  0.328320      0   1.00000 1.00636 0.083418
## 2  0.118576      1   0.67168 0.68689 0.056595
## 3  0.045810      2   0.55310 0.61092 0.049569
## 4  0.041761      3   0.50729 0.59117 0.047620
## 5  0.030918      4   0.46553 0.56251 0.044079
## 6  0.030731      5   0.43461 0.55924 0.045740
## 7  0.023733      6   0.40388 0.53676 0.044968
## 8  0.021595      7   0.38015 0.53695 0.043152
## 9  0.019798      8   0.35856 0.53188 0.043446
## 10 0.018870      9   0.33876 0.53251 0.042926
## 11 0.015480     10   0.31989 0.51882 0.042220
## 12 0.015420     11   0.30441 0.51052 0.041593
## 13 0.013580     12   0.28899 0.50808 0.041005
## 14 0.012335     13   0.27541 0.50372 0.041579
## 15 0.010000     14   0.26307 0.50122 0.041944
plotcp(carseats.tree)

min <- which.min(carseats.tree$cptable[,"xerror"])
print(min)
## 15 
## 15
cpt <- carseats.tree$cptable[which.min(carseats.tree$cptable[,"xerror"]),"CP"]
cpt
## [1] 0.01

con la validación cruzada obtenemos un tamaño de arbol igual a 13 y un cp de 0.01357961

ptree <- prune(carseats.tree,cp = cpt)
rpart.plot(ptree,box.palette = "RdBu",nn=TRUE)

pred1 <- predict(ptree,newdata = carseats.test)
m2 <- mean((pred1-carseats.test$Sales)^2)
m2
## [1] 5.463956

Se disminuye el MSE despues de que se poda.

(d) Utilice el enfoque de embolsado para analizar estos datos. ¿Qué prueba MSE obtienes? Use la función importance() para determinar qué variables son más importantes.

bag.carseats <- randomForest(Sales~.,carseats.train,mtry=(ncol(Carseats)-1) ,importance=TRUE)
pred2 <- predict(bag.carseats,newdata = carseats.test)
mean((pred2- carseats.test$Sales)^2)
## [1] 3.361658

notablemente mejora el MSE a 3.37627.

importance(bag.carseats)
##                 %IncMSE IncNodePurity
## CompPrice   26.67871283     183.54183
## Income       6.23358916      98.64888
## Advertising 13.57834499     110.41853
## Population   0.32338488      75.52763
## Price       68.47503499     599.01499
## ShelveLoc   77.26581305     725.54524
## Age         17.48598332     158.80999
## Education    0.01616168      57.99921
## Urban       -1.71716623      11.09547
## US           5.07615837      11.71181

Se puede concluir con lo anterior que Price,Shelveloc,CompPrice son variables importantes para este analisis.

(e) Utilice bosques aleatorios para analizar estos datos. ¿Qué prueba MSE obtienes? Use la función importance() para determinar qué variables son las más importantes. Describa el efecto de m, el número de variables consideradas en cada división, sobre la tasa de error obtenida.

bag.carseats <- randomForest(Sales ~ ., data = carseats.train, mtry = sqrt((ncol(Carseats) - 1)), importance = TRUE)
pred.bag <- predict(bag.carseats, newdata = carseats.test)
mean((pred.bag - carseats.test$Sales)^2)
## [1] 3.391034

Con lo anterior tenemos un MSE de 3.38

importance(bag.carseats)
##                %IncMSE IncNodePurity
## CompPrice   14.5403047     178.84227
## Income       2.6141067     154.02152
## Advertising 11.0107597     147.72102
## Population  -0.9294335     133.60346
## Price       40.3920222     478.45992
## ShelveLoc   50.2416174     559.76511
## Age         10.8775741     194.24705
## Education   -0.6566648      90.06653
## Urban       -0.4877354      18.95984
## US           3.6856553      23.97304

Con lo anterior cuncluimos que las mejores varibales predictoras para el analisis son las mismas que el numeral anterior.

###Punto 9

(a) Cree un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones, y un conjunto de prueba que contenga las observaciones restantes.

train <- sample(x = 1:nrow(OJ),size = 800)
OJ.train <- OJ[train,]
OJ.test <- OJ[-train,]

(b) Ajuste un árbol a los datos de entrenamiento, con Compra como respuesta y las otras variables como predictores. Use la función summary() para generar estadísticas resumidas sobre el árbol y describa los resultados obtenidos. ¿Cuál es la tasa de error de entrenamiento? ¿Cuántos nodos terminales tiene el árbol?

set.seed(0511)
library(tree)
tree.oj <- tree(Purchase~.,data = OJ.train)
summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH"        "SalePriceMM"    "SpecialCH"      "WeekofPurchase"
## [5] "ListPriceDiff"  "PctDiscMM"     
## Number of terminal nodes:  10 
## Residual mean deviance:  0.7071 = 558.6 / 790 
## Misclassification error rate: 0.1525 = 122 / 800

Este arbalo tiene 8 nodos terminales y una tasa de error de entrenamiento de 0.1612

(c) Escriba el nombre del objeto del árbol para obtener una salida de texto detallada. Elija uno de los nodos terminales e interprete la información que se muestra.

tree.oj
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 800 1063.000 CH ( 0.61875 0.38125 )  
##    2) LoyalCH < 0.5036 351  424.900 MM ( 0.29345 0.70655 )  
##      4) LoyalCH < 0.0356415 47    0.000 MM ( 0.00000 1.00000 ) *
##      5) LoyalCH > 0.0356415 304  389.300 MM ( 0.33882 0.66118 )  
##       10) SalePriceMM < 2.04 167  176.600 MM ( 0.22156 0.77844 )  
##         20) SpecialCH < 0.5 137  124.000 MM ( 0.16788 0.83212 ) *
##         21) SpecialCH > 0.5 30   41.460 MM ( 0.46667 0.53333 ) *
##       11) SalePriceMM > 2.04 137  189.700 MM ( 0.48175 0.51825 )  
##         22) LoyalCH < 0.168934 27   22.650 MM ( 0.14815 0.85185 ) *
##         23) LoyalCH > 0.168934 110  150.700 CH ( 0.56364 0.43636 )  
##           46) WeekofPurchase < 244.5 17   15.840 MM ( 0.17647 0.82353 ) *
##           47) WeekofPurchase > 244.5 93  122.100 CH ( 0.63441 0.36559 ) *
##    3) LoyalCH > 0.5036 449  341.700 CH ( 0.87305 0.12695 )  
##      6) LoyalCH < 0.764572 193  220.800 CH ( 0.74093 0.25907 )  
##       12) ListPriceDiff < 0.235 74  102.600 MM ( 0.50000 0.50000 )  
##         24) PctDiscMM < 0.196196 59   78.900 CH ( 0.61017 0.38983 ) *
##         25) PctDiscMM > 0.196196 15    7.348 MM ( 0.06667 0.93333 ) *
##       13) ListPriceDiff > 0.235 119   82.090 CH ( 0.89076 0.10924 ) *
##      7) LoyalCH > 0.764572 256   64.200 CH ( 0.97266 0.02734 ) *

Escogemos el nodo terminal 15, debido a que es un nodo terminal debido al asterisco.apreciamos que su criterio de división es LoyalCH > 0.764572 esta tiene 254 observaciones en este nodo y una desviación estandar 64.090. La predicción general para la rama CH,también podemos ver que el 97.24% de las observaciones toman el valor de CH, mientras que el 0.02756 tomal el valor de MM.

(d) Cree un diagrama del árbol e interprete los resultados.

plot(tree.oj)
text(tree.oj, pretty = 0)

Podemos apreciar que el nodo principal es LoyalCH si es menor o igual a 0.50 solo usa a su misma variables , pero si es menor utiliza a PriceDiff.

(e) Predecir la respuesta en los datos de prueba, y producir una matriz de confusión comparando las etiquetas de prueba con las etiquetas de prueba predichas. ¿Cuál es la tasa de error de prueba?

library(caret)
tree.pred <- predict(tree.oj, OJ.test, type = "class")

matrix<-confusionMatrix(tree.pred, OJ.test$Purchase)
matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  CH  MM
##         CH 138  33
##         MM  20  79
##                                           
##                Accuracy : 0.8037          
##                  95% CI : (0.7512, 0.8494)
##     No Information Rate : 0.5852          
##     P-Value [Acc > NIR] : 1.923e-14       
##                                           
##                   Kappa : 0.5887          
##                                           
##  Mcnemar's Test P-Value : 0.09929         
##                                           
##             Sensitivity : 0.8734          
##             Specificity : 0.7054          
##          Pos Pred Value : 0.8070          
##          Neg Pred Value : 0.7980          
##              Prevalence : 0.5852          
##          Detection Rate : 0.5111          
##    Detection Prevalence : 0.6333          
##       Balanced Accuracy : 0.7894          
##                                           
##        'Positive' Class : CH              
##

Observamos que la precisión del modelo es de 82.89%, siendo mas preciso para predecir a CH con un 90.65% mientras que para predecir MM sólo es de un 70.91%. (f) Aplique la función cv.tree() al conjunto de entrenamiento para determinar el tamaño óptimo del árbol.

cv.oj <- cv.tree(tree.oj,FUN = prune.misclass)
cv.oj
## $size
## [1] 10  9  6  2  1
## 
## $dev
## [1] 142 144 151 157 305
## 
## $k
## [1]       -Inf   0.000000   4.333333   6.250000 145.000000
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

(g) Produzca un diagrama con el tamaño del árbol en el eje xy la tasa de error de clasificación con validación cruzada en el eje y.

library(ggplot2)
library(reshape2)
valores <- data.frame(Size=cv.oj$size,Dev=cv.oj$dev)
theme_set(theme_bw())
ggplot(data = valores, aes(x=Size, y=Dev)) + 
  geom_line(color="red", linetype = "dashed") +
  geom_point() +
  scale_x_discrete(limits=c(1:9))

(h) ¿Qué tamaño de árbol corresponde a la tasa de error de clasificación con validación cruzada más baja? 5 es el tamaño con menor error

(i) Produzca un árbol podado correspondiente al tamaño óptimo del árbol obtenido mediante validación cruzada. Si la validación cruzada no conduce a la selección de un árbol podado, cree un árbol podado con cinco nodos terminales.

prune.oj <- prune.misclass(tree.oj,best = 5)
plot(prune.oj)
text(prune.oj,pretty = 0)

(j) Compare las tasas de error de entrenamiento entre los árboles podados y no podados. ¿Cuál es más alto?

summary(prune.oj)
## 
## Classification tree:
## snip.tree(tree = tree.oj, nodes = c(10L, 3L))
## Variables actually used in tree construction:
## [1] "LoyalCH"        "SalePriceMM"    "WeekofPurchase"
## Number of terminal nodes:  6 
## Residual mean deviance:  0.8552 = 679 / 794 
## Misclassification error rate: 0.1688 = 135 / 800
summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH"        "SalePriceMM"    "SpecialCH"      "WeekofPurchase"
## [5] "ListPriceDiff"  "PctDiscMM"     
## Number of terminal nodes:  10 
## Residual mean deviance:  0.7071 = 558.6 / 790 
## Misclassification error rate: 0.1525 = 122 / 800

observamos que hay una diferencia en el erro mayor el arbol podado. (k) Compare las tasas de error de prueba entre los árboles podados y no podados. ¿Cuál es más alto?

prune.pred <- predict(prune.oj, OJ.test, type = "class")

matrix<-confusionMatrix(prune.pred, OJ.test$Purchase)
matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  CH  MM
##         CH 140  39
##         MM  18  73
##                                          
##                Accuracy : 0.7889         
##                  95% CI : (0.7353, 0.836)
##     No Information Rate : 0.5852         
##     P-Value [Acc > NIR] : 1.161e-12      
##                                          
##                   Kappa : 0.553          
##                                          
##  Mcnemar's Test P-Value : 0.008071       
##                                          
##             Sensitivity : 0.8861         
##             Specificity : 0.6518         
##          Pos Pred Value : 0.7821         
##          Neg Pred Value : 0.8022         
##              Prevalence : 0.5852         
##          Detection Rate : 0.5185         
##    Detection Prevalence : 0.6630         
##       Balanced Accuracy : 0.7689         
##                                          
##        'Positive' Class : CH             
##

observamos que bajo la presicion del modelo notablemente con respecto al árbol sin podar.

###Punto 10

(a) Elimine las observaciones para las que se desconoce la información salarial, y luego transforme logarítmicamente los salarios.

Hitters <- na.omit(Hitters)
Hitters$Salary <- log(Hitters$Salary)

(b) Cree un conjunto de entrenamiento que consista en las primeras 200 observaciones, y un conjunto de prueba que consista en las observaciones restantes.

train <- 1:200
Hitters.train <- Hitters[train,]
Hitters.test <- Hitters[-train,]

(c) Realice un refuerzo en el conjunto de entrenamiento con 1,000 árboles para un rango de valores del parámetro de contracción \(\lambda\). Produzca un gráfico con diferentes valores de contracción en el eje xy el conjunto de entrenamiento MSE correspondiente en el eje y.

library(gbm)
## Loaded gbm 2.1.5
set.seed(0511)
lambdas = seq(0.001, 0.3, 0.005)
mse <- rep(0, length(lambdas))
for (i in 1:length(lambdas)) {
  boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
    pred.train <- predict(boost.hitters, Hitters.train, n.trees = 1000)
    mse[i] <- mean((pred.train - Hitters.train$Salary)^2)
}
library(ggplot2)

mse_graph <- data.frame(lambdas)
mse_graph["MSE"]<- mse
names(mse_graph) <- c("Lambdas", "MSE")

theme_set(theme_bw()) 
ggplot(data = mse_graph, aes(x = Lambdas, y = MSE)) +
  geom_line(color="red", linetype = "dashed") +
  geom_point()

(d) Produzca un gráfico con diferentes valores de contracción en el eje xy el conjunto de prueba MSE correspondiente en el eje y.

mse2 <- rep(NA, length(lambdas))
for (i in 1:length(lambdas)) {
    boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
    yhat <- predict(boost.hitters, Hitters.test, n.trees = 1000)
    mse2[i] <- mean((yhat - Hitters.test$Salary)^2)
}
library(ggplot2)

mse_graph2 <- data.frame(lambdas)
mse_graph2["MSE"]<- mse2
names(mse_graph2) <- c("Lambdas", "MSE")

theme_set(theme_bw()) 
ggplot(data = mse_graph2, aes(x = Lambdas, y = MSE)) +
  geom_line(color="red", linetype = "dashed") +
  geom_point() 

min(mse2)
## [1] 0.2461295
lambdas[which.min(mse2)]
## [1] 0.131

(e) Compare la prueba MSE de refuerzo con la prueba MSE que resulta de aplicar dos de los enfoques de regresión vistos en los Capítulos 3 y 6.

mod8.10 <- lm(Salary~.,data = Hitters.train)
pred8.10 <- predict(mod8.10,Hitters.test)
mean((pred8.10-Hitters.test$Salary)^2)
## [1] 0.4917959
library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
## 
##     R2
## The following object is masked from 'package:stats':
## 
##     loadings
pcr.fit=pcr(Salary~., data=Hitters.train , scale=TRUE , validation ="CV")
validationplot(pcr.fit ,val.type="MSEP")

pcr.pred=predict (pcr.fit ,Hitters.test,ncomp =1)
mean((pcr.pred -Hitters.test$Salary)^2)
## [1] 0.4661183

observamos que el MSE de la regresion lineal y de las componentes principales es mayor el MSE por el metodo Boosting

(f) ¿Qué variables parecen ser los predictores más importantes en el modelo impulsado?

boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[which.min(mse2)])
summary(boost.hitters)

##                 var    rel.inf
## CAtBat       CAtBat 22.1891184
## PutOuts     PutOuts  8.1156374
## CRBI           CRBI  7.7992242
## CWalks       CWalks  7.2267797
## Walks         Walks  6.9703139
## Assists     Assists  6.0854321
## CHmRun       CHmRun  5.2576243
## Hits           Hits  5.2422222
## Years         Years  4.8380328
## CRuns         CRuns  4.6471959
## RBI             RBI  4.6430518
## AtBat         AtBat  3.5618987
## HmRun         HmRun  3.5243625
## Errors       Errors  3.0632824
## CHits         CHits  2.7467221
## Runs           Runs  2.7155911
## Division   Division  0.6010006
## NewLeague NewLeague  0.5563958
## League       League  0.2161140

En la tabla se observa que la variable mas importante es CatBat.

(g) Ahora aplique el embolsado al conjunto de entrenamiento. ¿Cuál es el conjunto de prueba MSE para este enfoque?

bag.hitters <- randomForest(Salary ~ ., data = Hitters.train, mtry = 19)
yhat.bag <- predict(bag.hitters, newdata = Hitters.test)
mean((yhat.bag - Hitters.test$Salary)^2)
## [1] 0.228878

El MSE es de 0.2339 es mucho mejor que los otros metodos por el momento. ###Punto 11

(a) Cree un conjunto de entrenamiento que consista en las primeras 1,000 observaciones,y un conjunto de prueba que consta de las observaciones restantes.

(b) Ajuste un modelo de refuerzo al conjunto de entrenamiento con Compra como respuesta y las otras variables como predictores. Use 1,000 árboles y un valor de contracción de 0.01. ¿Qué predictores parecen ser los más importantes?

(c) Utilice el modelo de refuerzo para predecir la respuesta en los datos de la prueba. Predecir que una persona realizará una compra si la probabilidad estimada de compra es superior al 20%. Forme una matriz de confusión. ¿Qué fracción de las personas que predijeron hacer una compra en realidad la hacen? ¿Cómo se compara esto con los resultados obtenidos al aplicar KNN o regresión logística a este conjunto de datos?

###Punto 12

Capitulo IX

###Punto 4

library(e1071)
set.seed(1)
x <- rnorm(100)
y <- 4 * x^3 + 1 + rnorm(100)
class <- sample(100, 50)
y[class] <- y[class] + 3
y[-class] <- y[-class] - 3
plot(x[class], y[class], col = "green", xlab = "X", ylab = "Y", ylim = c(-10, 25), xlim = c(-1.5, 2.0))
points(x[-class], y[-class], col = "blue")

Clasificador de soporte en los datos de entrenamiento

z <- rep(-1, 100)
z[class] <- 1
data <- data.frame(x = x, y = y, z = as.factor(z))
train <- sample(100, 50)
data.train <- data[train, ]
data.test <- data[-train, ]
svm.linear <- svm(z ~ ., data = data.train, kernel = "linear", cost = 10)
plot(svm.linear, data.train)

table(predict = predict(svm.linear, data.train), truth = data.train$z)
##        truth
## predict -1  1
##      -1 21  1
##      1   4 24

El clasificador de vectores de soporte comete 4 errores en los datos de entrenamiento.

Maquina de vectores de soporte con un núcleo polinomial.

svm.poly <- svm(z ~ ., data = data.train, kernel = "polynomial", cost = 10)
plot(svm.poly, data.train)

table(predict = predict(svm.poly, data.train), truth = data.train$z)
##        truth
## predict -1  1
##      -1 23  1
##      1   2 24

La máquina de vectores de soporte comete 2 errores en los datos de entrenamiento.

svm.radial <- svm(z ~ ., data = data.train, kernel = "radial", gamma = 1, cost = 10)
plot(svm.radial, data.train)

table(predict = predict(svm.radial, data.train), truth = data.train$z)
##        truth
## predict -1  1
##      -1 25  0
##      1   0 25

La máquina con un núcleo radial comete 0 errores en los datos de entrenamiento.

Aplicando a los datos de prueba.

plot(svm.linear, data.test)

table(predict = predict(svm.linear, data.test), truth = data.test$z)
##        truth
## predict -1  1
##      -1 24  0
##      1   1 25
plot(svm.poly, data.test)

table(predict = predict(svm.poly, data.test), truth = data.test$z)
##        truth
## predict -1  1
##      -1 23  1
##      1   2 24
plot(svm.radial, data.test)

table(predict = predict(svm.radial, data.test), truth = data.test$z)
##        truth
## predict -1  1
##      -1 19  0
##      1   6 25

Como vemos en los resultados anteriores, las máquinas de vectores de soporte lineal, polinomial y radial clasifican, respectivamente 4, 2 y 0 observaciones de forma incorrecta. Entonces, el núcleo radial es el mejor modelo en este caso.

###Punto 5

  • Hemos visto que podemos ajustar un SVM con un núcleo no lineal para realizar la clasificación utilizando un límite de decisión no lineal. Ahora veremos que también podemos obtener un límite de decisión no lineal al realizar una regresión logística utilizando transformaciones no lineales de las características.

(a) Genere un conjunto de datos con n = 500 y p = 2, de modo que las observaciones pertenezcan a dos clases con un límite de decisión cuadrático entre ellas. Por ejemplo, puede hacer esto de la siguiente manera:

> x1 = runif (500) -0.5

> x2 = runif (500) -0.5

> y = 1 * (x1 ^ 2-x2 ^ 2> 0)

(b) Grafique las observaciones, coloreadas de acuerdo con sus etiquetas de clase. Su diagrama debe mostrar X1 en el eje xy X2 en el eje y.

(c) Ajuste un modelo de regresión logística a los datos, utilizando X1 y X2 como predictores.

(d) Aplique este modelo a los datos de entrenamiento para obtener una etiqueta de clase predicha para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas. El límite de decisión debe ser lineal.

(e) Ahora ajuste un modelo de regresión logística a los datos utilizando funciones no lineales de X1 y X2 como predictores (por ejemplo, X12, X1 × X2, log (X2), etc.).

(f) Aplique este modelo a los datos de entrenamiento para obtener una etiqueta de clase predicha para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas. El límite de decisión debe ser obviamente no lineal. Si no es así, repita (a) - (e) hasta que encuentre un ejemplo en el que las etiquetas de clase predichas sean obviamente no lineales.

(g) Ajuste un clasificador de vector de soporte a los datos con X1 y X2 como predictores. Obtenga una predicción de clase para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas.

(h) Ajuste un SVM usando un núcleo no lineal a los datos. Obtenga una predicción de clase para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas.

(i) Comente sus resultados.

###Punto 6

  • Al final de la Sección 9.6.1, se afirma que, en el caso de los datos que son apenas linealmente separables, un clasificador de vectores de soporte con un pequeño valor de costo que clasifica erróneamente un par de observaciones de entrenamiento puede funcionar mejor en los datos de prueba que uno con un enorme valor de costo que no clasifica erróneamente ninguna observación de capacitación. Ahora investigará este reclamo.

(a) Genere datos de dos clases con p = 2 de tal manera que las clases sean apenas separables linealmente.

(b) Calcule las tasas de error de validación cruzada para los clasificadores de vectores de soporte con un rango de valores de costo. ¿Cuántos errores de capacitación están mal clasificados para cada valor de costo considerado, y cómo se relaciona esto con los errores de validación cruzada obtenidos?

(c) Genere un conjunto de datos de prueba apropiado y calcule los errores de prueba correspondientes a cada uno de los valores de costo considerados. ¿Qué valor del costo conduce a la menor cantidad de errores de prueba, y cómo se compara esto con los valores de costo que producen la menor cantidad de errores de capacitación y la menor cantidad de errores de validación cruzada?

(d) Discuta sus resultados.

###Punto 7

  • En este problema, utilizará enfoques de vectores de soporte para predecir si un automóvil determinado obtiene un millaje de gasolina alto o bajo en función del conjunto de datos automático.

(a) Cree una variable binaria que tome un 1 para los automóviles con millaje de gasolina por encima de la mediana, y un 0 para automóviles con millaje de gasolina por debajo de la mediana.

(b) Ajuste un clasificador de vector de soporte a los datos con varios valores de costo, para predecir si un automóvil obtiene un millaje de gasolina alto o bajo. Informe los errores de validación cruzada asociados con diferentes valores de este parámetro. Comenta tus resultados.

(c) Ahora repita (b), esta vez usando SVM con núcleos de base radial y polinomial, con diferentes valores de gamma y grado y costo. Comenta tus resultados.

(d) Haga algunos gráficos para respaldar sus afirmaciones en (b) y (c). Sugerencia: En el laboratorio, utilizamos la función plot() para objetos svm solo en casos con p = 2. Cuando p> 2, puede usar la función plot() para crear gráficos que muestren pares de variables a la vez. Esencialmente, en lugar de escribir

> plot (svmfit, dat)

**donde svmfit contiene su modelo ajustado y dat es un marco de datos que contiene sus datos, puede escribir*

> plot (svmfit, dat, x1∼x4)

para graficar solo las variables primera y cuarta. Sin embargo, debe reemplazar x1 y x4 con los nombres correctos de las variables. Para obtener más información, escriba ?Plot.svm.

###Punto 8

  • Este problema involucra el conjunto de datos de OJ que es parte del paquete ISLR.

(a) Cree un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones, y un conjunto de prueba que contenga las observaciones restantes.

(b) Ajuste un clasificador de vector de soporte a los datos de entrenamiento usando cost = 0.01, con Compra como respuesta y las otras variables como predictores. Use la función summary() para generar estadísticas resumidas y describir los resultados obtenidos.

(c) ¿Cuáles son las tasas de error de capacitación y prueba?

(d) Use la función tune() para seleccionar un costo óptimo. Considere values en el rango de 0.01 a 10.

(e) Calcule las tasas de error de entrenamiento y prueba usando este nuevo valor por el cost

(f) Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo radial Use el valor predeterminado para gamma.

(g) Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo polinomial Establecer grado = 2.

(h) En general, ¿qué enfoque parece dar los mejores resultados con estos datos?